geom = readxl::read_xlsx("geometriePiece/ResultsAddedForExcel.xlsx",col_names=TRUE,skip=1)
New names:
• `` -> `...1`
colnames(geom)[1] = "type"
geom$type = as.factor(geom$type)
geom = split(geom,geom$type)$P225M909
Xgeom = geom[,2:4]
Y = geom[,5:ncol(geom)]
Gap = seq(30,7,by=-0.5)
matplot(Gap,t(Y),type="l",col = "steelblue",lwd=0.25,lty=1,xlab="Gap (mm)",ylab="Efforts (N)")
NAcurve = sort(unique(which(is.na(Y),arr.ind = TRUE)[,1]))
for (i in NAcurve){Y[i,which(is.na(Y[i,]))] = Y[i,min(which(is.na(Y[i,])))-1]}
# summary(Y) # OK, plus de donnees manquantes
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
gcv = 1:21
for (i in 1:21){
lambda = exp(i-10)
fdparTemp = fdPar(splbasis,Lfdobj = 2,lambda=lambda)
smoothdata = smooth.basis(argvals=Gap,t(Y),fdParobj = fdparTemp)
gcv[i] = mean(smoothdata$gcv)
}
lambda = exp(which.min(gcv)-10)
fdparTemp = fdPar(splbasis,Lfdobj = 2,lambda=lambda)
smoothdata = smooth.basis(argvals=Gap,t(Y),fdParobj = fdparTemp)
# plotfit.fd(t(Y),Gap,smoothdata$fd,pch=20,cex=0.5) #on est tres proche de l'interpolation
Ycoef = t(smoothdata$fd$coefs)
dim(Ycoef) # 14 coefficients pour n individus
[1] 498 14
#data = cbind(Ycoef,Xgeom)
#K=3;maxits=100;tol=1e-4; q = 3;p=14;nx=3
#est0 = estimates(cbind(Ycoef,Xgeom),K=3,par_init=NULL,maxits=300,
# tol=1e-4,
# q = 3,
# p=14,
# nx=3,
# cste = TRUE,
# verbose=TRUE)
#save(est0,file="geometriePiece/resultsGeom-31052023.Rdata")
load("geometriePiece/resultsGeom-31052023.Rdata")
Les résultats sont contenus dans l’objet est0
couleur = c("coral1","palegreen3","steelblue")
matplot(Gap,t(Y),type="l",col=couleur[est0$G],lwd=0.4,lty=1,xlab="Gap (mm)",ylab="Efforts (N)")
par(mfrow=c(1,3))
for (i in 1:length(est0$piik)){
matplot(Gap,t(Y[which(est0$G==i),]),type="l",col=couleur[i],ylim=c(0,42),lwd=0.4,lty=1,xlab="Gap (mm)",ylab="Efforts (N)")
}
On sélectionne d’abord les individus (courbes) d’un groupe particulier (tel que trouvé par MixPCA). On se concentre ici sur le groupe 1
grpe = 1 # a fixer suivant le groupe sur lequel on se concentre
index.gpe = which(est0$G==grpe)
print(paste("Le groupe ", grpe," contient ", length(index.gpe), " individus"))
[1] "Le groupe 1 contient 110 individus"
matplot(Gap,t(Y)[,index.gpe],type="l",col="steelblue",lty=1,lwd=0.4,xlab="Gap (mm)",ylab="Efforts (N)")
#reconstruction des courbes moyennes fonctionnelles
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
mufct = sapply(est0$mu,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
#ajout moyenne groupe 2
lines(Gap,mufct[,grpe],col="red",lwd=1)
On prend les courbes du groupe choisi, que l’on projette sur notre base de splines, ainsi que les géométries associées. On standardise les deux matrices (est-ce qu’il faut centrer uniquement?)
YcoefG = scale(Ycoef[index.gpe,])
XgeomG = scale(Xgeom[index.gpe,])
Les données de geometrie contiennent, dans l’ordre : * Thickness * Height * Shape
crosscor = matcor(YcoefG,XgeomG)
img.matcor(crosscor,type=2)
On se place dans un modèle du type
\[ Y = \mu_{Y} + Q_{Y} \alpha + \varepsilon \]
et
\[ X = \mu_{X} + Q_{X} \alpha + \eta \]
cc1 <- cc(XgeomG,YcoefG)
barplot(cc1$cor, xlab = "Dimension", ylab = "Canonical correlations", ylim = c(0,1))
cc1$cor
[1] 0.9987295 0.9793919 0.9360975
Les 3 corrélations canoniques sont plutot proches de 1
# Affichage des corrélations canoniques
cor(cc1$scores$xscores,cc1$scores$yscores)
[,1] [,2] [,3]
[1,] 9.987295e-01 -4.585389e-16 -5.155542e-16
[2,] 7.673963e-18 9.793919e-01 -6.345879e-17
[3,] 9.007057e-16 9.305821e-15 9.360975e-01
dim(cc1$scores$xscores)
[1] 110 3
Les scores selon X et selon Y (notés \(\alpha\)) sont très fortement corrélés et indépendant
La matrice \(Q_{X}\) est donnée par
cc1$xcoef
[,1] [,2] [,3]
Thickness -0.9364973 0.02020858 0.4268935
Height 0.1784473 1.03719849 -0.2303829
Shape 0.2744826 -0.18789032 1.0104856
t(cc1$xcoef)%*%cc1$xcoef
[,1] [,2] [,3]
[1,] 0.9842113 0.1145874 -0.1635351
[2,] 0.1145874 1.1114919 -0.4201863
[3,] -0.1635351 -0.4201863 1.2563954
cc1$xcoef%*%t(cc1$xcoef)
Thickness Height Shape
Thickness 1.0596736 -0.2445041 0.1705205
Height -0.2445041 1.1607004 -0.3786974
Shape 0.1705205 -0.3786974 1.1317245
xcoef1 = cc1$xcoef
Cette matrice nous renseigne sur la contribution des variables de géométrie à chaque axe factoriel.
La matrice \(Q_{Y}\) est donnée par
Qy = cc1$ycoef
Qy
[,1] [,2] [,3]
bspl4.1 0.04981210 -0.45107308 -0.7668100
bspl4.2 -0.05527291 -0.70640730 -2.3186012
bspl4.3 -0.41249359 -0.27723428 -2.4572998
bspl4.4 -0.26682343 0.99493213 -0.7637618
bspl4.5 -0.14870508 -0.66441240 6.8631132
bspl4.6 -0.15504243 -0.71888441 2.3721017
bspl4.7 0.05636064 1.86620349 -2.4302618
bspl4.8 -0.04081794 0.03082379 -1.4813699
bspl4.9 -0.08066014 -0.19333120 -0.3458932
bspl4.10 -0.06740741 0.19068015 0.2732237
bspl4.11 -0.03838069 0.40307386 0.3060224
bspl4.12 0.13905242 1.31726846 1.2849540
bspl4.13 0.29149592 2.42407301 2.6380758
bspl4.14 0.14690757 1.20261944 1.5097830
#t(cc1$ycoef)%*%cc1$ycoef
#V = t(cc1$scores$yscores)%*%YcoefG
#cor(t(V))
Elle correspond à nos composantes principales fonctionnelles. En “revenant” dans le domaine fonctionnel, on peut voir si celles-ci sont interprétables en tant que telles.
Qfct = apply(Qy,2,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
Qfct1 = Qfct
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l")
On représente (à vocation d’interprétation graphique)
\[ \mu_Y(t) \pm \alpha \hat{\psi}_j (t), \qquad j=1,\ldots,q\] où \(\alpha\) peut être choisi plus ou moins “arbitrairement” (normalement, on utilise la valeur propre associée pour moduler le coeff mais ici, ce n’est pas clair pour moi)
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,42),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qfct[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qfct[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
Il n’est pas évident de voir sur cette représentation si les composantes principales contenues dans \(Q_Y\) ont un sens fonctionnel.
Qvarimax = varimax(Qfct)
par(mfrow=c(1,2))
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l",main="sans VARIMAX")
matplot(Gap,Qvarimax$loadings,col=c("coral1","palegreen3","steelblue"),type="l",main="avec VARIMAX")
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,42),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qvarimax$loadings[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qvarimax$loadings[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
La différence n’est pas complètement claire…
Par ailleurs, on est limité à 3 composantes principales étant donné qu’on a 3 covariables de géométrie. Philosophiquement, ce n’est sans doute pas étonnant puisqu’on a au final des courbes qui n’ont que 3 degrés de liberté.
On sélectionne d’abord les individus (courbes) d’un groupe particulier (tel que trouvé par MixPCA). On se concentre ici sur le groupe 2
grpe = 2 # a fixer suivant le groupe sur lequel on se concentre
index.gpe = which(est0$G==grpe)
print(paste("Le groupe ", grpe," contient ", length(index.gpe), " individus"))
[1] "Le groupe 2 contient 179 individus"
matplot(Gap,t(Y)[,index.gpe],type="l",col="steelblue",lty=1,lwd=0.4,xlab="Gap (mm)",ylab="Efforts (N)")
#reconstruction des courbes moyennes fonctionnelles
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
mufct = sapply(est0$mu,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
#ajout moyenne groupe 2
lines(Gap,mufct[,grpe],col="red",lwd=1)
On prend les courbes du groupe choisi, que l’on projette sur notre base de splines, ainsi que les géométries associées. On standardise les deux matrices (est-ce qu’il faut centrer uniquement?)
YcoefG = scale(Ycoef[index.gpe,])
XgeomG = scale(Xgeom[index.gpe,])
Les données de geometrie contiennent, dans l’ordre : * Thickness * Height * Shape
crosscor = matcor(YcoefG[,1:9],XgeomG)
img.matcor(crosscor,type=2)
On se place dans un modèle du type
\[ Y = \mu_{Y} + Q_{Y} \alpha + \varepsilon \]
et
\[ X = \mu_{X} + Q_{X} \alpha + \eta \]
cc1 <- cc(XgeomG,YcoefG[,1:9]) # on enlève les 5 derniers coefs (singulrite matrice)
barplot(cc1$cor, xlab = "Dimension", ylab = "Canonical correlations", ylim = c(0,1))
cc1$cor
[1] 0.9870316 0.8924068 0.4784289
Les 3 corrélations canoniques sont moins proches de 1 ici (interpretation à reflechir)
# Affichage des corrélations canoniques
cor(cc1$scores$xscores,cc1$scores$yscores)
[,1] [,2] [,3]
[1,] 9.870316e-01 1.377484e-16 -1.320109e-15
[2,] -4.273553e-17 8.924068e-01 1.458884e-16
[3,] -1.461115e-16 1.758456e-17 4.784289e-01
dim(cc1$scores$xscores)
[1] 179 3
Les scores selon X et selon Y (notés \(\alpha\)) sont très fortement corrélés et indépendant
La matrice \(Q_{X}\) est donnée par
cc1$xcoef
[,1] [,2] [,3]
Thickness 0.8410331 -0.7386838 0.3001075
Height 0.4434320 0.9675255 -0.3651595
Shape -0.4890760 0.3397731 0.9320679
#t(cc1$xcoef)%*%cc1$xcoef # orthogonalite?
#cc1$xcoef%*%t(cc1$xcoef) # orthogonalite?
xcoef2 = cc1$xcoef
Cette matrice nous renseigne sur la contribution des variables de géométrie à chaque axe factoriel. Elle semble a priori assez différente de celle du groupe 1 meme si on retrouve un deuxième “axe” lié à la hauteur et un troisième “axe” lié à shape. Note : elle n’est toujours pas orthogonale
La matrice \(Q_{Y}\) est donnée par
Qy = cc1$ycoef
Qy
[,1] [,2] [,3]
bspl4.1 0.08938044 -0.52163186 1.4863663
bspl4.2 0.06420464 -0.21975696 -1.7278458
bspl4.3 0.06382005 0.45138160 0.8170389
bspl4.4 0.13458351 -2.38261643 2.3933461
bspl4.5 0.79450371 0.52431956 -2.9886906
bspl4.6 -0.12820481 2.11949191 -0.8735578
bspl4.7 -0.07424761 0.21352924 0.8470028
bspl4.8 0.08737024 -0.08318324 0.3348507
bspl4.9 0.05717563 -0.25571302 -0.1567053
#t(cc1$ycoef)%*%cc1$ycoef
#V = t(cc1$scores$yscores)%*%YcoefG
#cor(t(V))
Elle correspond à nos composantes principales fonctionnelles. En “revenant” dans le domaine fonctionnel, on peut voir si celles-ci sont interprétables en tant que telles.
Ycoeffin = matrix(apply(YcoefG,2,mean)[10:14],nrow=5,ncol=3) # on reprend la moyenne des coefs pour la fin des courbes (grands gaps)
Qfct = apply(rbind(Qy,Ycoeffin),2,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
Qfct2 = Qfct
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l")
On représente (à vocation d’interprétation graphique)
\[ \mu_Y(t) \pm \alpha \hat{\psi}_j (t), \qquad j=1,\ldots,q\] où \(\alpha\) peut être choisi plus ou moins “arbitrairement” (normalement, on utilise la valeur propre associée pour moduler le coeff mais ici, ce n’est pas clair pour moi)
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,42),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qfct[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qfct[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
On peut interpréter les courbes en regard de la matrice
xcoef. Ces courbes nous renseigne sur les parties impactées
par la géométrie en fonction de l’axe factoriel considéré. Par exemple,
le deuxième axe de xcoef est positivement lié à
Height mais négativement à Thickness. La
reconstruction fonctionnelle du deuxième axe ycoef est
cohérente avec cela car on observe que les fortes valeurs positives
(grande hauteur/faible épaisseur) montre une croissance de la courbe
plus précoce que pour les fortes valeurs négatives (petites
hauteurs/forte épaisseur)
On sélectionne d’abord les individus (courbes) d’un groupe particulier (tel que trouvé par MixPCA). On se concentre ici sur le groupe 3
grpe = 3 # a fixer suivant le groupe sur lequel on se concentre
index.gpe = which(est0$G==grpe)
print(paste("Le groupe ", grpe," contient ", length(index.gpe), " individus"))
[1] "Le groupe 3 contient 209 individus"
matplot(Gap,t(Y)[,index.gpe],type="l",col="steelblue",lty=1,lwd=0.4,xlab="Gap (mm)",ylab="Efforts (N)")
#reconstruction des courbes moyennes fonctionnelles
splbasis = create.bspline.basis(rangeval=range(Gap),norder=4,breaks=seq(min(Gap),max(Gap),length=12))
mufct = sapply(est0$mu,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
#ajout moyenne groupe 2
lines(Gap,mufct[,grpe],col="red",lwd=1)
On prend les courbes du groupe choisi, que l’on projette sur notre base de splines, ainsi que les géométries associées. On voit ici que le groupe 3 n’est pas très homogène, il contient des courbes avec des comportements très différents. De même que pour les deux autres groupes, on standardise les deux matrices (est-ce qu’il faut centrer uniquement?)
YcoefG = scale(Ycoef[index.gpe,])
XgeomG = scale(Xgeom[index.gpe,])
Les données de geometrie contiennent, dans l’ordre : * Thickness * Height * Shape
crosscor = matcor(YcoefG,XgeomG)
img.matcor(crosscor,type=2)
Les corrélations croisées présentent un comportement relativement différents dans ce groupe (en 3 blocs au lieu de 2), faisant echo à l’hétérogénéité des courbes dans le groupe.
On se place dans un modèle du type
\[ Y = \mu_{Y} + Q_{Y} \alpha + \varepsilon \]
et
\[ X = \mu_{X} + Q_{X} \alpha + \eta \]
cc1 <- cc(XgeomG,YcoefG)
barplot(cc1$cor, xlab = "Dimension", ylab = "Canonical correlations", ylim = c(0,1))
cc1$cor
[1] 0.9760231 0.9679074 0.3567894
Pas de souci de singularité dans ce groupe. La 3eme corrélation canonique est relativement faible ici (interpretation à reflechir)
# Affichage des corrélations canoniques
cor(cc1$scores$xscores,cc1$scores$yscores)
[,1] [,2] [,3]
[1,] 9.760231e-01 -1.909368e-16 -5.523116e-16
[2,] 6.306779e-18 9.679074e-01 4.844452e-16
[3,] 1.436371e-16 3.311742e-17 3.567894e-01
#dim(cc1$scores$xscores)
Les scores selon X et selon Y (notés \(\alpha\)) sont très fortement corrélés et indépendant (sauf la 3eme)
La matrice \(Q_{X}\) est donnée par
cc1$xcoef
[,1] [,2] [,3]
Thickness 0.6909726 0.6703926 0.4054039
Height -0.4893503 0.9174937 0.0650506
Shape -0.2755391 -0.1476526 0.9551033
#t(cc1$xcoef)%*%cc1$xcoef # orthogonalite?
#cc1$xcoef%*%t(cc1$xcoef) # orthogonalite?
xcoef3 = cc1$xcoef
Cette matrice nous renseigne sur la contribution des variables de géométrie à chaque axe factoriel. On retrouve une matrice légèrement différente des autres groupes même si l’axe 2 est toujours lié positivement à la hauteur, l’axe 1 est lié à l’épaisseur positivement et négativement à la hauteur. L’axe 3 reste lié à shape.
La matrice \(Q_{Y}\) est donnée par
Qy = cc1$ycoef
Qy
[,1] [,2] [,3]
bspl4.1 0.94397278 -0.150832261 -1.147420132
bspl4.2 -0.42956322 -0.807084948 -0.783172061
bspl4.3 0.22854101 0.414421492 1.654343278
bspl4.4 0.57389907 0.753557408 0.498022864
bspl4.5 -0.46776895 0.157110102 0.009861507
bspl4.6 -0.51906239 0.042890446 2.158460635
bspl4.7 -0.26390244 0.033380575 -0.464832741
bspl4.8 0.63205519 0.269904302 -4.024685529
bspl4.9 -0.05048885 0.284290003 1.711892056
bspl4.10 -0.10579046 0.008809529 0.471992567
bspl4.11 -0.03258642 0.109354682 0.313024771
bspl4.12 -0.01188870 -0.086751591 2.634823744
bspl4.13 0.03706895 -0.210887782 3.993030181
bspl4.14 0.04946376 -0.093288446 1.255825534
#t(cc1$ycoef)%*%cc1$ycoef
#V = t(cc1$scores$yscores)%*%YcoefG
#cor(t(V))
Elle correspond à nos composantes principales fonctionnelles. En “revenant” dans le domaine fonctionnel, on peut voir si celles-ci sont interprétables en tant que telles.
Qfct = apply(Qy,2,FUN=function(x) {ft =fd(x,splbasis) ; return(eval.fd(Gap,ft))})
Qfct3 = Qfct
matplot(Gap,Qfct,col=c("coral1","palegreen3","steelblue"),type="l")
On représente (à vocation d’interprétation graphique)
\[ \mu_Y(t) \pm \alpha \hat{\psi}_j (t), \qquad j=1,\ldots,q\] où \(\alpha\) peut être choisi plus ou moins “arbitrairement” (normalement, on utilise la valeur propre associée pour moduler le coeff mais ici, ce n’est pas clair pour moi)
for(i in 1:3){
matplot(Gap,t(Y)[,index.gpe],type="l",col="gray",ylab="Efforts",ylim=c(0,20),lwd=0.5)
lines(Gap,mufct[,grpe],col=1)
lines(Gap,mufct[,grpe]+3*Qfct[,i],type="p",pch="+",col=couleur[i])
lines(Gap,mufct[,grpe]-3*Qfct[,i],type="p",pch="-",col=couleur[i],lwd=2)
}
Le même type d’interprétation peut prévaloir ici mais les résultats sont plus difficiles à interpréter, le groupe 3 semblant rassembler des courbes très diverses (groupe un peu fourre-tout)
par(mfrow=c(1,3))
image.plot(matrix(rev(as.vector(t(xcoef1))),3,3)[3:1,], main="Groupe 1",xaxt="n",yaxt="n")
axis(side=2,at=c(0,0.5,1),labels=c("Shape","Height","Thickness"),tick=F)
axis(side=1,at=c(0,0.5,1),labels=c("Axe1","Axe2","Axe3"),tick=F)
image.plot(matrix(rev(as.vector(t(xcoef2))),3,3)[3:1,], main="Groupe 2",xaxt="n",yaxt="n")
axis(side=2,at=c(0,0.5,1),labels=c("Shape","Height","Thickness"),tick=F)
axis(side=1,at=c(0,0.5,1),labels=c("Axe1","Axe2","Axe3"),tick=F)
image.plot(matrix(rev(as.vector(t(xcoef3))),3,3)[3:1,], main="Groupe 3",xaxt="n",yaxt="n")
axis(side=2,at=c(0,0.5,1),labels=c("Shape","Height","Thickness"),tick=F)
axis(side=1,at=c(0,0.5,1),labels=c("Axe1","Axe2","Axe3"),tick=F)
On a des matrices relativement différentes dans les 3 groupes, même si on les regarde à transposition près. NB : dans le groupe 2, la CCA est appliquée en enlevant les 5 derniers coefficients splines des courbes.
Au niveau des axes principaux “fonctionnels”, on a aussi des axes assez différents
par(mfrow=c(1,3))
matplot(Gap,Qfct1,col=c("coral1","palegreen3","steelblue"),type="l")
matplot(Gap,Qfct2,col=c("coral1","palegreen3","steelblue"),type="l")
matplot(Gap,Qfct3,col=c("coral1","palegreen3","steelblue"),type="l")
Il est plus difficile d’expliquer les différences entre les groupes dans cette analyse double. Mais ici, on est parti avec des groupes définis à l’avance, il faudra voir ce que cela donne si on construit des groupes en se basant sur la CCA (ou proche)